<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
   "http://www.w3.org/TR/html4/strict.dtd">

<html>
<head>
  <title></title>
  <meta http-equiv="content-type" content="text/html; charset=latin1">
  <style type="text/css">
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
pre { line-height: 125%; }
body .hll { background-color: #ffffcc }
body  { background: #f8f8f8; }
body .c { color: #008800; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #AA22FF; font-weight: bold } /* Keyword */
body .o { color: #666666 } /* Operator */
body .cm { color: #008800; font-style: italic } /* Comment.Multiline */
body .cp { color: #008800 } /* Comment.Preproc */
body .c1 { color: #008800; font-style: italic } /* Comment.Single */
body .cs { color: #008800; font-weight: bold } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #888888 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0044DD } /* Generic.Traceback */
body .kc { color: #AA22FF; font-weight: bold } /* Keyword.Constant */
body .kd { color: #AA22FF; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #AA22FF; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #AA22FF } /* Keyword.Pseudo */
body .kr { color: #AA22FF; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #00BB00; font-weight: bold } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #BB4444 } /* Literal.String */
body .na { color: #BB4444 } /* Name.Attribute */
body .nb { color: #AA22FF } /* Name.Builtin */
body .nc { color: #0000FF } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #00A000 } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #008000; font-weight: bold } /* Name.Tag */
body .nv { color: #B8860B } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sb { color: #BB4444 } /* Literal.String.Backtick */
body .sc { color: #BB4444 } /* Literal.String.Char */
body .sd { color: #BB4444; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #BB4444 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #BB4444 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #008000 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #BB4444 } /* Literal.String.Single */
body .ss { color: #B8860B } /* Literal.String.Symbol */
body .bp { color: #AA22FF } /* Name.Builtin.Pseudo */
body .vc { color: #B8860B } /* Name.Variable.Class */
body .vg { color: #B8860B } /* Name.Variable.Global */
body .vi { color: #B8860B } /* Name.Variable.Instance */
body .il { color: #666666 } /* Literal.Number.Integer.Long */

  </style>
</head>
<body>
<h2></h2>

<div class="highlight"><pre><span class="c">#from libmyatomism import *</span>
<span class="c">#from math import *</span>
<span class="c">#export LD_PRELOAD=/usr/lib/libmpi.so python</span>
<span class="kn">import</span> <span class="nn">array</span> 
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">plt</span>
<span class="kn">from</span> <span class="nn">libatomism</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">math</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">mpl_toolkits.mplot3d</span> <span class="kn">import</span> <span class="n">Axes3D</span>
<span class="kn">from</span> <span class="nn">matplotlib.mlab</span> <span class="kn">import</span> <span class="n">griddata</span>
<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">cm</span>


<span class="sd">&quot;&quot;&quot; this function plot the 2D PES isocontours, for your </span>
<span class="sd">first lecture you should go directly to the main body</span>
<span class="sd">of the code.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="k">def</span> <span class="nf">plotIsoPesForApproach</span><span class="p">(</span> <span class="n">kinop</span><span class="p">,</span> <span class="n">repro</span><span class="p">,</span> <span class="n">fig</span> <span class="p">,</span> <span class="n">nfig</span><span class="p">):</span>

   <span class="n">zmat</span> <span class="o">=</span> <span class="n">kinop</span><span class="o">.</span><span class="n">getSystem</span><span class="p">()</span>
   <span class="n">gcoors</span> <span class="o">=</span> <span class="n">kinop</span><span class="o">.</span><span class="n">getCoordinates</span><span class="p">()</span>
   <span class="n">q0</span> <span class="o">=</span> <span class="n">gcoors</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
   <span class="n">q1</span> <span class="o">=</span> <span class="n">gcoors</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
   <span class="n">q0</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mf">2.</span><span class="p">,</span><span class="mf">2.</span><span class="p">,</span><span class="mf">5.</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">0.1</span><span class="p">)</span>
   <span class="n">q1</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">90.1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>

   <span class="sd">&quot;&quot;&quot; A linear sampler: scan from min to max like a &#39;for&#39; loop</span>
<span class="sd">   imbriqued for each activated coordinates (if Dq&gt;0)</span>
<span class="sd">   &quot;&quot;&quot;</span>
   <span class="n">scan</span><span class="o">=</span><span class="n">msSamplerLinear</span><span class="o">.</span><span class="n">New</span><span class="p">()</span>
  
   <span class="n">gcoors</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">.</span><span class="n">Dq</span><span class="o">=</span><span class="mi">0</span>
   <span class="n">gcoors</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span><span class="o">.</span><span class="n">Dq</span><span class="o">=</span><span class="mi">0</span>
   <span class="sd">&quot;&quot;&quot; the iso levels we want to plot &quot;&quot;&quot;</span>
   <span class="n">levels</span><span class="o">=</span><span class="p">[</span><span class="o">-</span><span class="mf">0.9</span><span class="p">,</span><span class="o">-</span><span class="mf">0.75</span><span class="p">,</span><span class="o">-</span><span class="mf">0.5</span><span class="p">,</span><span class="o">-</span><span class="mf">0.25</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mf">2.5</span><span class="p">,</span><span class="mi">5</span><span class="p">,</span><span class="mi">10</span><span class="p">,</span><span class="mi">20</span><span class="p">,</span><span class="mi">40</span><span class="p">,</span><span class="mi">50</span><span class="p">,</span><span class="mi">60</span><span class="p">,</span><span class="mi">70</span><span class="p">,</span><span class="mi">80</span><span class="p">]</span>
   <span class="n">zeroindice</span><span class="o">=</span><span class="mi">4</span>
   
   <span class="sd">&quot;&quot;&quot; create the array which will contain the values of the 2D PES</span>
<span class="sd">   &quot;&quot;&quot;</span>
   <span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="n">q0</span><span class="o">.</span><span class="n">Min</span><span class="p">,</span><span class="n">q0</span><span class="o">.</span><span class="n">Max</span><span class="p">,</span><span class="n">q0</span><span class="o">.</span><span class="n">Dq</span><span class="p">)</span>
   <span class="n">Y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="n">q1</span><span class="o">.</span><span class="n">Min</span><span class="p">,</span><span class="n">q1</span><span class="o">.</span><span class="n">Max</span><span class="p">,</span><span class="n">q1</span><span class="o">.</span><span class="n">Dq</span><span class="p">)</span>
   <span class="n">Z</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">Y</span><span class="p">),</span><span class="nb">len</span><span class="p">(</span><span class="n">X</span><span class="p">)</span> <span class="p">))</span>  
   <span class="n">X</span><span class="p">,</span> <span class="n">Y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">meshgrid</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">)</span>

   <span class="n">scan</span><span class="o">.</span><span class="n">begin</span><span class="p">()</span>
   <span class="k">while</span><span class="p">(</span> <span class="n">scan</span><span class="o">.</span><span class="n">getNextPoint</span><span class="p">(</span> <span class="n">gcoors</span> <span class="p">)</span> <span class="p">):</span>
  
       <span class="sd">&quot;&quot;&quot; set the coordinates and the atoms&#39; position</span>
<span class="sd">       &quot;&quot;&quot;</span>    
       <span class="n">kinop</span><span class="o">.</span><span class="n">setDynamicDofs</span><span class="p">()</span>
       <span class="n">zmat</span><span class="o">.</span><span class="n">computeCartCoordinates</span><span class="p">()</span>    

       <span class="sd">&quot;&quot;&quot; Here we get the values of the qij_repro</span>
<span class="sd">       &quot;&quot;&quot;</span>
       <span class="n">vec</span> <span class="o">=</span> <span class="n">VectorOfDouble</span><span class="p">()</span>
       <span class="n">vec</span><span class="o">.</span><span class="n">extend</span><span class="p">([</span><span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">4</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">),</span><span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">5</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">),</span><span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">4</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">),</span><span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">5</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">)])</span>
       <span class="n">a</span> <span class="o">=</span>  <span class="n">repro</span><span class="o">.</span><span class="n">evaluate</span><span class="p">(</span><span class="n">vec</span><span class="p">)</span>

       <span class="k">if</span><span class="p">(</span> <span class="p">(</span><span class="n">a</span><span class="o">&gt;</span><span class="n">Emax</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">a</span><span class="o">&lt;-</span><span class="mi">10</span><span class="p">)</span> <span class="p">):</span>
           <span class="n">a</span><span class="o">=</span><span class="n">Emax</span>
       <span class="n">Z</span><span class="p">[</span><span class="n">q1</span><span class="o">.</span><span class="n">indice</span><span class="p">()][</span><span class="n">q0</span><span class="o">.</span><span class="n">indice</span><span class="p">()]</span> <span class="o">=</span> <span class="n">a</span>
   <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">   plot</span>
<span class="sd">   &quot;&quot;&quot;</span>
   <span class="n">plt</span><span class="o">=</span><span class="n">fig</span><span class="o">.</span><span class="n">add_subplot</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">nfig</span><span class="p">)</span>
   <span class="n">CS</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">contour</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">,</span> <span class="n">Z</span><span class="p">,</span><span class="n">levels</span><span class="p">,</span> <span class="n">colors</span><span class="o">=</span><span class="s">&#39;k&#39;</span><span class="p">)</span>
   <span class="n">zc</span> <span class="o">=</span> <span class="n">CS</span><span class="o">.</span><span class="n">collections</span><span class="p">[</span><span class="n">zeroindice</span><span class="p">]</span> 
   <span class="n">zc</span><span class="o">.</span><span class="n">set_linewidth</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>
   <span class="n">plt</span><span class="o">.</span><span class="n">clabel</span><span class="p">(</span><span class="n">CS</span><span class="p">,</span> <span class="n">inline</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="mi">12</span><span class="p">)</span>
   <span class="n">plt</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;q0 [Angstrom]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
   <span class="n">plt</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;q1 [Degree]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
   <span class="n">plt</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>


<span class="sd">&quot;&quot;&quot; main body of the code</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">msLogger</span><span class="o">.</span><span class="n">setPriorityMin</span><span class="p">(</span><span class="n">Priority</span><span class="o">.</span><span class="n">ERROR</span><span class="p">)</span>

<span class="n">unit</span><span class="o">=</span><span class="n">msUnitsManager</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom Degree amu kcal/mol ps&quot;</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Declaration of the z-matrix</span>
<span class="sd">X </span>
<span class="sd">N 0 0.55 </span>
<span class="sd">N 0 0.55        1  179.9  </span>
<span class="sd">X 0 5.   (=q0)  1  90.   (=q1)   2   0. </span>
<span class="sd">N 3 0.55        0  90.   (=q2)   1   180. (=q3)</span>
<span class="sd">N 3 0.55        4  179.9         0   180.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">zmat</span> <span class="o">=</span> <span class="n">msZMat</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;zmat&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="s">&quot;X N 0 0.55 N 0 0.55 1 179.9 X 0 5. 1 90. 2 0. N 3 0.55 0 90. 1 180. N 3 0.55 4 179.9 0 180.&quot;</span><span class="p">)</span>

<span class="n">system</span><span class="o">=</span><span class="n">msSystem</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">addEntity</span><span class="p">(</span><span class="n">zmat</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Declaration of the dynamic variables of the z-matrix, </span>
<span class="sd">those who will control the relative orientation of the two N2.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">q0</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">5</span><span class="p">,</span><span class="mf">0.5</span><span class="p">,</span><span class="mi">6</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">0.1</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q0&quot;</span><span class="p">)</span>
<span class="n">q1</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Degree&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">120</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">179.9</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q1&quot;</span><span class="p">)</span>
<span class="n">q2</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Degree&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">120</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">179.9</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q2&quot;</span><span class="p">)</span>
<span class="n">q3</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Degree&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">120</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">360</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q3&quot;</span><span class="p">)</span>

<span class="n">gcoors</span> <span class="o">=</span> <span class="n">msGeneralizedCoordinates</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q0</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q1</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q2</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q3</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">The kinetic operator: it does the link between generalized </span>
<span class="sd">coordinates and geometry</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">kinop</span><span class="o">=</span><span class="n">msKineticOperator</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="n">gcoors</span><span class="p">,</span><span class="n">system</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">the following is defining the kinetic operator with dynamic DoF of type f(q)=q;</span>
<span class="sd">the first definition is the same as doing, for instance:</span>
<span class="sd">  kinfct0=msParser.New(&quot;Angstrom&quot;).setCoordinates(gcoors).setExpression(&quot;q0&quot;).setId(&quot;f=q0&quot;)</span>
<span class="sd">  kinop.addDynamicDof(zmat,3,kinfct0)</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">kinop</span><span class="o">.</span><span class="n">addDynamicDofLinear</span><span class="p">(</span><span class="n">zmat</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="n">q0</span><span class="p">)</span><span class="o">.</span><span class="n">addDynamicDofLinear</span><span class="p">(</span><span class="n">zmat</span><span class="p">,</span><span class="mi">4</span><span class="p">,</span><span class="n">q1</span><span class="p">)</span><span class="o">.</span><span class="n">addDynamicDofLinear</span><span class="p">(</span><span class="n">zmat</span><span class="p">,</span><span class="mi">7</span><span class="p">,</span><span class="n">q2</span><span class="p">)</span><span class="o">.</span><span class="n">addDynamicDofLinear</span><span class="p">(</span><span class="n">zmat</span><span class="p">,</span><span class="mi">8</span><span class="p">,</span><span class="n">q3</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; </span>
<span class="sd">Load the database ( format of each line:  q0 q1 q2 q3 q4 E_in_hartree E_in_kcal/mol ) </span>
<span class="sd">we want to fit E_in_kcal/mol (0 of energy is 2*E(N2) )</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="k">with</span> <span class="nb">open</span><span class="p">(</span><span class="s">&#39;Database.txt&#39;</span><span class="p">)</span> <span class="k">as</span> <span class="n">f</span><span class="p">:</span>
    <span class="n">w</span> <span class="o">=</span> <span class="n">f</span><span class="o">.</span><span class="n">readline</span><span class="p">()</span><span class="o">.</span><span class="n">split</span><span class="p">()</span> <span class="c"># read first line</span>
    <span class="n">array</span> <span class="o">=</span> <span class="p">[]</span>
    <span class="k">for</span> <span class="n">line</span> <span class="ow">in</span> <span class="n">f</span><span class="p">:</span>           <span class="c"># read rest of lines</span>
        <span class="n">array</span><span class="o">.</span><span class="n">append</span><span class="p">([</span><span class="nb">float</span><span class="p">(</span><span class="n">x</span><span class="p">)</span> <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="n">line</span><span class="o">.</span><span class="n">split</span><span class="p">()])</span>

<span class="sd">&quot;&quot;&quot; </span>
<span class="sd">Define the variable of the fitting function: </span>
<span class="sd">this is not the gcoors because we want to use variable of type 1/(R*R)</span>
<span class="sd">where R is the atoms separation. For instance, q02 refer to the separation</span>
<span class="sd"> between the nitrogen atom number O and number 2.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">q02_repro</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom^-2&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mf">1.e6</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">0.1</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q0_repro&quot;</span><span class="p">)</span>
<span class="n">q03_repro</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom^-2&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mf">1.e6</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">0.1</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q1_repro&quot;</span><span class="p">)</span>
<span class="n">q12_repro</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom^-2&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mf">1.e6</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">0.1</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q2_repro&quot;</span><span class="p">)</span>
<span class="n">q13_repro</span> <span class="o">=</span> <span class="n">msScalarVariable</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;Angstrom^-2&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mf">1.e6</span><span class="p">,</span><span class="mf">0.1</span><span class="p">,</span><span class="mf">0.1</span><span class="p">)</span><span class="o">.</span><span class="n">setId</span><span class="p">(</span><span class="s">&quot;q3_repro&quot;</span><span class="p">)</span>
<span class="n">gcoors_repro</span> <span class="o">=</span> <span class="n">msGeneralizedCoordinates</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="n">unit</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q02_repro</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q03_repro</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q12_repro</span><span class="p">)</span><span class="o">.</span><span class="n">addVar</span><span class="p">(</span><span class="n">q13_repro</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Repromodel is a type of scalar function, we initialize it </span>
<span class="sd">with the unit of the returned value, and set the coordinates after that.</span>
<span class="sd">&quot;&quot;&quot;</span> 
<span class="n">repro</span> <span class="o">=</span><span class="n">msRepromodel</span><span class="o">.</span><span class="n">New</span><span class="p">(</span><span class="s">&quot;kcal/mol&quot;</span><span class="p">)</span><span class="o">.</span><span class="n">setCoordinates</span><span class="p">(</span><span class="n">gcoors_repro</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">This is the order of the development, and the maximum energy accounted</span>
<span class="sd">&quot;&quot;&quot;</span> 
<span class="n">N</span> <span class="o">=</span> <span class="mi">5</span>
<span class="n">Emax</span><span class="o">=</span><span class="mi">80</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Here we add the data to fit into the repromodel.</span>
<span class="sd">The database inside the repromodel is a scalar function of </span>
<span class="sd">type dirac.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="k">for</span> <span class="n">values</span> <span class="ow">in</span> <span class="n">array</span><span class="p">:</span>

    <span class="n">vec</span> <span class="o">=</span> <span class="n">VectorOfDouble</span><span class="p">()</span>
    <span class="n">vec</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span> <span class="n">arg</span> <span class="k">for</span> <span class="n">arg</span> <span class="ow">in</span> <span class="n">values</span><span class="p">)</span>
    <span class="sd">&quot;&quot;&quot; set the coordinates and the atoms&#39; position</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">gcoors</span><span class="o">.</span><span class="n">setCoordinates</span><span class="p">(</span><span class="n">vec</span><span class="p">)</span>
    <span class="n">kinop</span><span class="o">.</span><span class="n">setDynamicDofs</span><span class="p">()</span>
    <span class="n">system</span><span class="o">.</span><span class="n">computeCartCoordinates</span><span class="p">()</span>
    <span class="sd">&quot;&quot;&quot; Here we get the values of the qij_repro</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">vec</span> <span class="o">=</span> <span class="n">VectorOfDouble</span><span class="p">()</span>
    <span class="n">vec</span><span class="o">.</span><span class="n">extend</span><span class="p">([</span> <span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">4</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">)</span> <span class="p">,</span> <span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">5</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">)</span> <span class="p">,</span> <span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">4</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">)</span> <span class="p">,</span> <span class="n">zmat</span><span class="o">.</span><span class="n">separation</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span><span class="mi">5</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="p">)])</span>
    <span class="k">if</span><span class="p">(</span> <span class="n">values</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span><span class="o">&lt;</span><span class="n">Emax</span> <span class="p">):</span>
        <span class="n">repro</span><span class="o">.</span><span class="n">addData</span><span class="p">(</span> <span class="n">vec</span><span class="p">,</span> <span class="n">values</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span> <span class="p">)</span>
        
<span class="n">msLogger</span><span class="o">.</span><span class="n">setPriorityMin</span><span class="p">(</span><span class="n">Priority</span><span class="o">.</span><span class="n">INFO</span><span class="p">)</span>
<span class="n">msLogger</span><span class="o">.</span><span class="n">setVerboseLevel</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; Create the orthonormal basis set of order N</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">repro</span><span class="o">.</span><span class="n">setBasisSet</span><span class="p">(</span><span class="n">N</span><span class="p">)</span>
<span class="n">msLogger</span><span class="o">.</span><span class="n">setPriorityMin</span><span class="p">(</span><span class="n">Priority</span><span class="o">.</span><span class="n">ERROR</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; fit the function: projection of the database function (i.e dirac function)</span>
<span class="sd">on the basis set. </span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">repro</span><span class="o">.</span><span class="n">fit</span><span class="p">()</span>

<span class="sd">&quot;&quot;&quot; After this point, this is only for plotting the results ...</span>

<span class="sd">return a sampler containing the scenario </span>
<span class="sd">of the database providen</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">sampler</span> <span class="o">=</span> <span class="n">repro</span><span class="o">.</span><span class="n">getSampler</span><span class="p">()</span>

<span class="sd">&quot;&quot;&quot; get the database: it is a scalar</span>
<span class="sd">function of type dirac, each time evaluated </span>
<span class="sd">for the sampler current position. see the </span>
<span class="sd">tutorial page for details.</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">field</span> <span class="o">=</span> <span class="n">repro</span><span class="o">.</span><span class="n">getField</span><span class="p">()</span>

<span class="n">y</span><span class="o">=</span><span class="p">[]</span>
<span class="n">x</span><span class="o">=</span><span class="p">[]</span>
<span class="n">n</span><span class="o">=</span><span class="mi">0</span>
<span class="n">sigma</span><span class="o">=</span><span class="mi">0</span>

<span class="sd">&quot;&quot;&quot; begin the sampling &quot;&quot;&quot;</span>
<span class="n">sampler</span><span class="o">.</span><span class="n">begin</span><span class="p">()</span>

<span class="sd">&quot;&quot;&quot;comparison between extact and fitted values</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="k">while</span><span class="p">(</span><span class="n">sampler</span><span class="o">.</span><span class="n">getNextPoint</span><span class="p">(</span><span class="n">gcoors_repro</span><span class="p">)):</span>
     <span class="n">x</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">field</span><span class="o">.</span><span class="n">evaluate</span><span class="p">())</span>
     <span class="n">y</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">repro</span><span class="o">.</span><span class="n">evaluate</span><span class="p">())</span>
     <span class="n">sigma</span><span class="o">+=</span><span class="nb">abs</span><span class="p">(</span><span class="n">field</span><span class="o">.</span><span class="n">evaluate</span><span class="p">()</span> <span class="o">-</span> <span class="n">repro</span><span class="o">.</span><span class="n">evaluate</span><span class="p">())</span>
     <span class="n">n</span><span class="o">=</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span>
<span class="k">print</span> <span class="s">&quot;sigma=&quot;</span><span class="p">,</span><span class="n">sigma</span><span class="o">/</span><span class="n">n</span>


<span class="n">fig</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">plt</span><span class="o">=</span><span class="n">fig</span><span class="o">.</span><span class="n">add_subplot</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">axis</span><span class="p">((</span><span class="o">-</span><span class="mi">5</span><span class="p">,</span><span class="n">Emax</span><span class="p">,</span><span class="o">-</span><span class="mi">5</span><span class="p">,</span><span class="n">Emax</span><span class="p">))</span>
<span class="n">plt</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s">&quot;Exact energy [kcal/mol]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s">&quot;Fitted energy [kcal/mol]&quot;</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">tick_params</span><span class="p">(</span><span class="n">labelsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>

<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="s">&quot;+&quot;</span><span class="p">)</span>


<span class="sd">&quot;&quot;&quot; First approach</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">q2</span><span class="o">.</span><span class="n">Value</span><span class="o">=</span><span class="mf">179.9</span>
<span class="n">q3</span><span class="o">.</span><span class="n">Value</span><span class="o">=</span><span class="mf">0.</span>
<span class="n">plotIsoPesForApproach</span><span class="p">(</span> <span class="n">kinop</span><span class="p">,</span> <span class="n">repro</span><span class="p">,</span> <span class="n">fig</span> <span class="p">,</span> <span class="mi">2</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; same for the second approach</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">q2</span><span class="o">.</span><span class="n">Value</span><span class="o">=</span><span class="mf">90.</span>
<span class="n">plotIsoPesForApproach</span><span class="p">(</span> <span class="n">kinop</span><span class="p">,</span> <span class="n">repro</span><span class="p">,</span> <span class="n">fig</span> <span class="p">,</span> <span class="mi">3</span><span class="p">)</span>

<span class="sd">&quot;&quot;&quot; and the third approach</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="n">q2</span><span class="o">.</span><span class="n">Value</span><span class="o">=</span><span class="mf">90.</span>
<span class="n">q3</span><span class="o">.</span><span class="n">Value</span><span class="o">=</span><span class="mf">90.</span>

<span class="n">plotIsoPesForApproach</span><span class="p">(</span> <span class="n">kinop</span><span class="p">,</span> <span class="n">repro</span><span class="p">,</span> <span class="n">fig</span> <span class="p">,</span> <span class="mi">4</span><span class="p">)</span>

<span class="n">fig</span><span class="o">.</span><span class="n">set_size_inches</span><span class="p">(</span><span class="mi">20</span><span class="p">,</span><span class="mi">20</span><span class="p">)</span>
<span class="n">fig</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="s">&#39;N4_0.png&#39;</span><span class="p">)</span>

<span class="n">fig</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>  
</pre></div>
</body>
</html>
